 SUBROUTINE TRANSPOSE_XY_YZ(F,nxpsi,nxpse,nypsi,nypse,nzpsi,nzpse, &
     &                     Q,imin,imax,jmin,jmax,kmin,kmax,nx,ny,nz, &
     &                     comm3d,nproc,ndim,dims,coords, &
     &                     nxpi,nxpe,nypi,nype,nzpi,nzpe,ksize)

!This subroutine transposes the XY array decomposition into
!YZ decomposition

 !use MPI
 use MPI_VARIABLES, only : &
 & transpose_xy_yz_indic
 use TIMING_MOD, only : &
 & proc_time

 implicit none

 include 'mpif.h'


!Input variables 
 real(8), intent(in) :: F(nxpsi:nxpse,nypsi:nypse,nzpsi:nzpse)
 integer, intent(in) :: nxpsi,nxpse,nypsi,nypse,nzpsi,nzpse
 integer, intent(in) :: imin, imax
 integer, intent(in) :: jmin, jmax
 integer, intent(in) :: kmin, kmax
 integer, intent(in) :: nx, ny, nz 
 integer, intent(in) :: ndim, nproc
 integer, intent(in) :: comm3d
 integer, intent(in) :: coords(1:3)
 integer, intent(in) :: dims(1:3)
 integer, intent(in) :: nxpi(0:nproc-1)
 integer, intent(in) :: nxpe(0:nproc-1)
 integer, intent(in) :: nypi(0:nproc-1)
 integer, intent(in) :: nype(0:nproc-1)
 integer, intent(in) :: nzpi(0:nproc-1)
 integer, intent(in) :: nzpe(0:nproc-1)
 integer, intent(in) :: ksize(0:nproc-1)

!Output variables
 real(8), intent(out) :: Q(imin:imax,jmin:jmax,kmin:kmax)

!Local variables
 integer, allocatable, save :: maxbuf(:) ! The maximal size of send buffer for each process 

 integer :: i, j, k
 integer :: ki, icoord
 integer :: coordn(1:3)
 integer :: coordrecv(1:3)
 integer :: ibuf,kbuf,tag,rrank,rankrecv,ierr
 integer :: status(MPI_STATUS_SIZE,2*nproc)
 integer :: status1(MPI_STATUS_SIZE)
 integer, allocatable, save :: req2(:)
 integer :: myrank
 integer, save :: max_maxbuf

 real(8), allocatable, save :: sarr(:,:) ! Send array
 real(8), allocatable, save :: rarr(:,:) ! Receive arrayi
 real(8) :: time00 ! Work variable
 real(8) :: pr_time, pro_time

 logical, save :: firstcall = .true.

! pr_time = MPI_WTIME()
 call cpu_time(pr_time)

 call MPI_COMM_RANK(comm3d, myrank, ierr)

 if (firstcall) then
   allocate (maxbuf(0:nproc-1) )
   allocate (req2(1:2*nproc) )
   call MAXBUFCALC()
   max_maxbuf = 1
   do rrank = 0, nproc - 1
     max_maxbuf = max(max_maxbuf, maxbuf(rrank))
   enddo
   allocate (sarr(1:maxbuf(myrank), 0:nproc-1) )
   allocate (rarr(1:max_maxbuf,     0:nproc-1) )
 endif


 tag = 1
 coordn(2) = coords(2); coordrecv(2) = coords(2)
 coordn(3) = coords(3); coordrecv(3) = coords(3)

 ki = nzpi(myrank)
 if (coords(1) /= 0) then
   do icoord = 0, coords(1) - 1
     ki = ki + ksize(icoord)
   enddo
 endif

 coordrecv(1) = coords(1)

! time00 = MPI_WTIME()

 do icoord = coords(1), dims(1) - 1
   coordn(1) = icoord   
   call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
   call MPI_CART_RANK(COMM3D,coordrecv,rankrecv,ierr)
   ibuf = 1
   do k = ki, ki + ksize(icoord) - 1 ! Sending for coordn(1) = icoord process the ksize(icoord) portion 
     do j = nypi(myrank), nype(myrank)
       do i = nxpi(myrank), nxpe(myrank)
         sarr(ibuf,rrank) = F(i,j,k)
         ibuf = ibuf + 1
       enddo
     enddo
   enddo


!   if (firstcall) then
!   call MPI_send_init(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
!     &                rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

!   call MPI_recv_init(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
!     &                rankrecv,tag,comm3d,req2(coordrecv(1)+1),ierr)

!   endif

!   call MPI_START(req2(icoord+nproc+1), ierr)
!   call MPI_START(req2(coordrecv(1)+1), ierr)


   call MPI_ISEND(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
     &                rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

   call MPI_IRECV(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
     &                rankrecv,tag,comm3d,req2(coordrecv(1)+1),ierr)


   ki = ki + ksize(icoord)
   coordrecv(1) = coordrecv(1) - 1
   if (coordrecv(1) == -1) coordrecv(1) = dims(1) - 1
 enddo


 if (coords(1) /= 0) then
 ki = nzpi(myrank)
 do icoord = 0, coords(1) - 1
   coordn(1) = icoord
   call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
   call MPI_CART_RANK(COMM3D,coordrecv,rankrecv,ierr)
   ibuf = 1
   do k = ki, ki + ksize(icoord) - 1 ! Sending for coordn(1) = icoord process the ksize(icoord) portion 
     do j = nypi(myrank), nype(myrank)
       do i = nxpi(myrank), nxpe(myrank)
         sarr(ibuf,rrank) = F(i,j,k)
         ibuf = ibuf + 1
       enddo
     enddo
   enddo

!   if (firstcall) then
!   call MPI_send_init(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
!     &                rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

!   call MPI_recv_init(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
!     &                rankrecv,tag,comm3d,req2(coordrecv(1)+1),ierr)

!   endif

!   call MPI_START(req2(icoord+nproc+1), ierr)
!   call MPI_START(req2(coordrecv(1)+1), ierr)


   call MPI_ISEND(sarr(1,rrank),maxbuf(myrank),MPI_DOUBLE_PRECISION, &
     &                rrank,tag,comm3d,req2(icoord+nproc+1),ierr)

   call MPI_IRECV(rarr(1,rankrecv),maxbuf(rankrecv),MPI_DOUBLE_PRECISION, &
     &                rankrecv,tag,comm3d,req2(coordrecv(1)+1),ierr)


   ki = ki + ksize(icoord)
   coordrecv(1) = coordrecv(1) - 1
   if (coordrecv(1) == -1) coordrecv(1) = dims(1) - 1
 enddo
 endif

! time_(1) = time_(1) + MPI_WTIME() - time00

! time00 = MPI_WTIME()

 call MPI_WAITALL(dims(1), req2(1), status(1,1), ierr)

! time_(2) = time_(2) + MPI_WTIME() - time00

! time00 = MPI_WTIME()

 coordn(2) = coords(2)
 do icoord = 0, dims(1) - 1
   coordn(1) = icoord
   call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
   ibuf = 1
   do k = kmin, kmin + ksize(coords(1)) - 1
     do j = nypi(rrank), nype(rrank)
       do i = nxpi(rrank), nxpe(rrank)
         Q(i,j,k) = rarr(ibuf,rrank)
         ibuf = ibuf + 1
       enddo
     enddo
   enddo
 enddo

! time_(3) = time_(3) + MPI_WTIME() - time00

! time00 = MPI_WTIME()
 
 call MPI_WAITALL(dims(1), req2(nproc+1), status(1,1), ierr)

! time_(4) = time_(4) + MPI_WTIME() - time00

 call MPI_BARRIER(comm3d, ierr)

! proc_time(transpose_xy_yz_indic,2) = proc_time(transpose_xy_yz_indic,2) + &
! & MPI_WTIME() - pr_time

 call cpu_time(pro_time)
 proc_time(transpose_xy_yz_indic,2) = proc_time(transpose_xy_yz_indic,2) + &
 & pro_time - pr_time

 if (firstcall) firstcall = .false.
 RETURN
 contains


 SUBROUTINE MAXBUFCALC

 implicit none

 k = 1
 do icoord = 0, dims(1) - 1
   k = max(k, ksize(icoord))
 enddo

 do j = 0, dims(2) - 1
   do i = 0, dims(1) - 1
     coordn(1) = i
     coordn(2) = j
     coordn(3) = 0
     call MPI_CART_RANK(COMM3D,coordn,rrank,ierr)
     maxbuf(rrank) = (nxpe(rrank) - nxpi(rrank) + 1) * &
     &               (nype(rrank) - nypi(rrank) + 1) * k
   enddo
 enddo

 END SUBROUTINE MAXBUFCALC

 END SUBROUTINE TRANSPOSE_XY_YZ
